#----Basics----

rm(list = ls())

# Set working directory
setwd("~/Library/CloudStorage/Dropbox/Whitten_Kagalwala/TS/JOP/Acceptance/KW_Replication/GECM_DGP")

# Load required libraries
library(tidyverse)
library(haven)

# Load data
data <- read_dta('gecm.dta')
data$time <- 1:7

# Reshape data
phi1_bias <- data %>%
  select(c(time, starts_with("phi_bias"), starts_with("phi1_bias"))) %>%
  pivot_longer(phi_bias_gecm:phi1_bias_adl33, names_to = "model", values_to = "phi1_bias") %>%
  mutate(model = gsub("^ph.*_bias_(.*)", "\\1", model))
phi1_type1 <- data %>%
  select(c(time,matches("^phi_(.*)_t1$"), matches("^phi1_(.*)_t1$"))) %>%
  pivot_longer(phi_gecm_t1:phi1_adl33_t1, names_to = "model", values_to = "phi1_type1") %>%
  mutate(model = gsub("^ph.*_(.*)_t1", "\\1", model))
phi1_power <- data %>%
  select(c(time,matches("^phi_(.*)_power$"), matches("^phi1_(.*)_power$"))) %>%
  pivot_longer(phi_gecm_power:phi1_adl33_power, names_to = "model", values_to = "phi1_power") %>%
  mutate(model = gsub("^ph.*_(.*)_power", "\\1", model))
phi2_bias <- data %>%
  select(c(time, starts_with("phi2_bias"))) %>%
  pivot_longer(phi2_bias_adl22:phi2_bias_adl33, names_to = "model", values_to = "phi2_bias") %>%
  mutate(model = gsub("^ph.*_bias_(.*)", "\\1", model))
phi2_type1 <- data %>%
  select(c(time,matches("^phi2_(.*)_t1$"))) %>%
  pivot_longer(phi2_adl22_t1:phi2_adl33_t1, names_to = "model", values_to = "phi2_type1") %>%
  mutate(model = gsub("^ph.*_(.*)_t1", "\\1", model))
phi3_bias <- data %>%
  select(c(time, starts_with("phi3_bias"))) %>%
  pivot_longer(phi3_bias_adl33, names_to = "model", values_to = "phi3_bias") %>%
  mutate(model = gsub("^ph.*_bias_(.*)", "\\1", model))
phi3_type1 <- data %>%
  select(c(time,matches("^phi3_(.*)_t1$"))) %>%
  pivot_longer(phi3_adl33_t1, names_to = "model", values_to = "phi3_type1") %>%
  mutate(model = gsub("^ph.*_(.*)_t1", "\\1", model))
b1_bias <- data %>%
  select(c(time, starts_with("b1_bias"))) %>%
  pivot_longer(b1_bias_gecm:b1_bias_adl33, names_to = "model", values_to = "b1_bias") %>%
  mutate(model = gsub("^b1_bias_(.*)", "\\1", model))
b1_type1 <- data %>%
  select(c(time,matches("^b1_(.*)_t1$"))) %>%
  pivot_longer(b1_gecm_t1:b1_adl33_t1, names_to = "model", values_to = "b1_type1") %>%
  mutate(model = gsub("^b1_(.*)_t1", "\\1", model))
b1_power <- data %>%
  select(c(time,matches("^b1_(.*)_power$"))) %>%
  pivot_longer(b1_gecm_power:b1_adl33_power, names_to = "model", values_to = "b1_power") %>%
  mutate(model = gsub("^b1_(.*)_power", "\\1", model))
b2_bias <- data %>%
  select(c(time, starts_with("b2_bias"))) %>%
  pivot_longer(b2_bias_gecm:b2_bias_adl33, names_to = "model", values_to = "b2_bias") %>%
  mutate(model = gsub("^b2_bias_(.*)", "\\1", model))
b2_type1 <- data %>%
  select(c(time,matches("^b2_(.*)_t1$"))) %>%
  pivot_longer(b2_gecm_t1:b2_adl33_t1, names_to = "model", values_to = "b2_type1") %>%
  mutate(model = gsub("^b2_(.*)_t1", "\\1", model))
b2_power <- data %>%
  select(c(time,matches("^b2_(.*)_power$"))) %>%
  pivot_longer(b2_gecm_power:b2_adl33_power, names_to = "model", values_to = "b2_power") %>%
  mutate(model = gsub("^b2_(.*)_power", "\\1", model))
b3_bias <- data %>%
  select(c(time, starts_with("b3_bias"))) %>%
  pivot_longer(b3_bias_adl22:b3_bias_adl33, names_to = "model", values_to = "b3_bias") %>%
  mutate(model = gsub("^b3_bias_(.*)", "\\1", model))
b3_type1 <- data %>%
  select(c(time,matches("^b3_(.*)_t1$"))) %>%
  pivot_longer(b3_adl22_t1:b3_adl33_t1, names_to = "model", values_to = "b3_type1") %>%
  mutate(model = gsub("^b3_(.*)_t1", "\\1", model))
b4_bias <- data %>%
  select(c(time, starts_with("b4_bias"))) %>%
  pivot_longer(b4_bias_adl33, names_to = "model", values_to = "b4_bias") %>%
  mutate(model = gsub("^b4_bias_(.*)", "\\1", model))
b4_type1 <- data %>%
  select(c(time,matches("^b4_(.*)_t1$"))) %>%
  pivot_longer(b4_adl33_t1, names_to = "model", values_to = "b4_type1") %>%
  mutate(model = gsub("^b4_(.*)_t1", "\\1", model))
lrm_bias <- data %>%
  select(c(time, starts_with("lrm_bias"))) %>%
  pivot_longer(lrm_bias_gecm:lrm_bias_adl33, names_to = "model", values_to = "lrm_bias") %>%
  mutate(model = gsub("^lrm_bias_(.*)", "\\1", model))
lrm_type1 <- data %>%
  select(c(time,matches("^lrm_(.*)_t1$"))) %>%
  pivot_longer(lrm_gecm_t1:lrm_adl33_t1, names_to = "model", values_to = "lrm_type1") %>%
  mutate(model = gsub("^lrm_(.*)_t1", "\\1", model))
lrm_power <- data %>%
  select(c(time,matches("^lrm_(.*)_power$"))) %>%
  pivot_longer(lrm_gecm_power:lrm_adl33_power, names_to = "model", values_to = "lrm_power") %>%
  mutate(model = gsub("^lrm_(.*)_power", "\\1", model))
stationary_y <- data %>%
  select(c(time,stat_y)) 
ftest_y_t1 <- data %>%
  dplyr::select(c(time,matches("^ftest_y_(.*)_t1$"))) %>%
  pivot_longer(ftest_y_adl33_t1:ftest_y_adl22_t1, names_to = "model", values_to = "ftest_y_t1") %>%
  mutate(model = gsub("^ftest.*_(.*)_t1", "\\1", model))
ftest_x_t1 <- data %>%
  dplyr::select(c(time,matches("^ftest_x_(.*)_t1$"))) %>%
  pivot_longer(ftest_x_adl33_t1:ftest_x_adl22_t1, names_to = "model", values_to = "ftest_x_t1") %>%
  mutate(model = gsub("^ftest.*_(.*)_t1", "\\1", model))
ftest_y_power <- data %>%
  dplyr::select(c(time,matches("^ftest_y_(.*)_power$"))) %>%
  pivot_longer(ftest_y_adl33_power:ftest_y_adl22_power, names_to = "model", values_to = "ftest_y_power") %>%
  mutate(model = gsub("^ftest.*_(.*)_power", "\\1", model))
ftest_x_power <- data %>%
  dplyr::select(c(time,matches("^ftest_x_(.*)_power$"))) %>%
  pivot_longer(ftest_x_adl33_power:ftest_x_adl22_power, names_to = "model", values_to = "ftest_x_power") %>%
  mutate(model = gsub("^ftest.*_(.*)_power", "\\1", model))

#----y stationarity----
ggplot(stationary_y) +
  geom_hline(yintercept = 0.05, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = stat_y), position = position_dodge(width = 0.65)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 2, byrow = TRUE)) +
  labs(y = "Type 1 Error Rate", x = "T", shape = "", title = expression("Stationarity of y"[t])) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10)) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000"))
ggsave("GECM_staty_type1.pdf")

#----Ly Bias----
phi1_bias$model <- factor(phi1_bias$model, levels = c("gecm", "adl", "adl22", "adl33"), labels = c("gecm", "adl", "adl22", "adl33"))
ggplot(phi1_bias, aes(shape = model)) +
  geom_hline(yintercept = 0, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = phi1_bias), position = position_dodge(width = 0.75)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Bias", x = "T", shape = "", title = expression("Bias in the coefficient of y"[t-1])) + ylim(-0.5, 0.5) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("gecm", "adl", "adl22", "adl33"),
                     values = c(3, 4, 5, 6),
                     labels = c("GECM", "ADL(1,1)", "ADL(2,2)", "ADL(3,3)")) +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1) 
ggsave("GECM_phi1_bias.pdf")

#----Ly Type 1----
phi1_type1$model <- factor(phi1_type1$model, levels = c("gecm", "adl", "adl22", "adl33"), labels = c("gecm", "adl", "adl22", "adl33"))
ggplot(phi1_type1, aes(shape = model)) +
  geom_hline(yintercept = 0.05, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = phi1_type1), position = position_dodge(width = 0.65)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Type 1 Error Rate", x = "T", shape = "", title = expression("Type 1 error rate in the coefficient of y"[t-1])) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("gecm", "adl", "adl22", "adl33"),
                     values = c(3, 4, 5, 6),
                     labels = c("GECM", "ADL(1,1)", "ADL(2,2)", "ADL(3,3)")) +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1) 
ggsave("GECM_phi1_type1.pdf")

#----Ly Coverage----
phi1_type1 %>%
  mutate(phi1_coverage = 1 - phi1_type1) %>%
ggplot(aes(shape = model)) +
  geom_hline(yintercept = 0.95, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = phi1_coverage), position = position_dodge(width = 0.65)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Coverage Probability", x = "T", shape = "", title = expression("Coverage probability of the coefficient of y"[t-1])) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("gecm", "adl", "adl22", "adl33"),
                     values = c(3, 4, 5, 6),
                     labels = c("GECM", "ADL(1,1)", "ADL(2,2)", "ADL(3,3)")) +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1) 
ggsave("GECM_phi1_coverage.pdf")

#----Ly Power----
phi1_power$model <- factor(phi1_power$model, levels = c("gecm", "adl", "adl22", "adl33"), labels = c("gecm", "adl", "adl22", "adl33"))
ggplot(phi1_power, aes(shape = model)) +
  geom_hline(yintercept = 1, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = phi1_power), position = position_dodge(width = 0.65)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Power", x = "T", shape = "", title = expression("Power in the coefficient of y"[t-1])) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) +  
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("gecm", "adl", "adl22", "adl33"),
                     values = c(3, 4, 5, 6),
                     labels = c("GECM", "ADL(1,1)", "ADL(2,2)", "ADL(3,3)")) +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1) 
ggsave("GECM_phi1_power.pdf")

#----L2y Bias----
phi2_bias$model <- factor(phi2_bias$model, levels = c("adl22", "adl33"), labels = c("adl22", "adl33"))
ggplot(phi2_bias, aes(shape = model)) +
  geom_hline(yintercept = 0, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = phi2_bias), position = position_dodge(width = 0.35)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Bias", x = "T", shape = "", title = expression("Bias in the coefficient of y"[t-2])) + ylim(-0.5, 0.5) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("adl22", "adl33"),
                     values = c(5, 6),
                     labels = c("ADL(2,2)", "ADL(3,3)")) +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1) 
ggsave("GECM_phi2_bias.pdf")

#----L2y Type 1----
phi2_type1$model <- factor(phi2_type1$model, levels = c("adl22", "adl33"), labels = c("adl22", "adl33"))
ggplot(phi2_type1, aes(shape = model)) +
  geom_hline(yintercept = 0.05, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = phi2_type1), position = position_dodge(width = 0.35)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Type 1 Error Rate", x = "T", shape = "", title = expression("Type 1 error rate in the coefficient of y"[t-2])) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) +  
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("adl22", "adl33"),
                     values = c(5, 6),
                     labels = c("ADL(2,2)", "ADL(3,3)")) +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1) 
ggsave("GECM_phi2_type1.pdf")

#----L2y Coverage----
phi2_type1 %>%
  mutate(phi2_coverage = 1 - phi2_type1) %>%
ggplot(aes(shape = model)) +
  geom_hline(yintercept = 0.95, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = phi2_coverage), position = position_dodge(width = 0.35)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Coverage Probability", x = "T", shape = "", title = expression("Coverage probability of the coefficient of y"[t-2])) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) +  
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("adl22", "adl33"),
                     values = c(5, 6),
                     labels = c("ADL(2,2)", "ADL(3,3)")) +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1) 
ggsave("GECM_phi2_coverage.pdf")

#----L3y Bias----
ggplot(phi3_bias, aes(shape = model)) +
  geom_hline(yintercept = 0, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = phi3_bias), position = position_dodge(width = 0.35)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Bias", x = "T", shape = "", title = expression("Bias in the coefficient of y"[t-3])) + ylim(-0.5, 0.5) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10)) +
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("adl33"),
                     values = c(6),
                     labels = c("ADL(3,3)"))
ggsave("GECM_phi3_bias.pdf")

#----L3y Type 1----
ggplot(phi3_type1, aes(shape = model)) +
  geom_hline(yintercept = 0.05, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = phi3_type1), position = position_dodge(width = 0.35)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Type 1 Error Rate", x = "T", shape = "", title = expression("Type 1 error rate in the coefficient of y"[t-3])) + ylim(0, 1) + 
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10)) +
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("adl33"),
                     values = c(6),
                     labels = c("ADL(3,3)"))
ggsave("GECM_phi3_type1.pdf")

#----L3y Coverage----
phi3_type1 %>%
  mutate(phi3_coverage = 1 - phi3_type1) %>%
ggplot(aes(shape = model)) +
  geom_hline(yintercept = 0.95, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = phi3_coverage), position = position_dodge(width = 0.35)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Coverage Probability", x = "T", shape = "", title = expression("Coverage probability of the coefficient of y"[t-3])) + ylim(0, 1) + 
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10)) +
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("adl33"),
                     values = c(6),
                     labels = c("ADL(3,3)"))
ggsave("GECM_phi3_coverage.pdf")

#----x Bias----
b1_bias$model <- factor(b1_bias$model, levels = c("gecm", "adl", "adl22", "adl33"), labels =  c("gecm", "adl", "adl22", "adl33"))
ggplot(b1_bias, aes(shape = model)) +
  geom_hline(yintercept = 0, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = b1_bias), position = position_dodge(width = 0.75)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Bias", x = "T", shape = "", title = expression("Bias in the coefficient of x"[t])) + ylim(-0.5, 0.5) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) +  
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("gecm", "adl", "adl22", "adl33"),
                     values = c(3, 4, 5, 6),
                     labels = c("GECM", "ADL(1,1)", "ADL(2,2)", "ADL(3,3)")) +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1) 
ggsave("GECM_b1_bias.pdf")

#----x Type 1----
b1_type1$model <- factor(b1_type1$model, levels = c("gecm", "adl", "adl22", "adl33"), labels =  c("gecm", "adl", "adl22", "adl33"))
ggplot(b1_type1, aes(shape = model)) +
  geom_hline(yintercept = 0.05, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = b1_type1), position = position_dodge(width = 0.75)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Type 1 Error Rate", x = "T", shape = "", title = expression("Type 1 error rate in the coefficient of x"[t])) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("gecm", "adl", "adl22", "adl33"),
                     values = c(3, 4, 5, 6),
                     labels = c("GECM", "ADL(1,1)", "ADL(2,2)", "ADL(3,3)")) +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1) 
ggsave("GECM_b1_type1.pdf")

#----x Coverage----
b1_type1 %>%
  mutate(b1_coverage = 1 - b1_type1) %>%
  ggplot(aes(shape = model)) +
  geom_hline(yintercept = 0.95, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = b1_coverage), position = position_dodge(width = 0.75)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Coverage Probability", x = "T", shape = "", title = expression("Coverage probability of the coefficient of x"[t])) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("gecm", "adl", "adl22", "adl33"),
                     values = c(3, 4, 5, 6),
                     labels = c("GECM", "ADL(1,1)", "ADL(2,2)", "ADL(3,3)")) +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1) 
ggsave("GECM_b1_coverage.pdf")

#----x Power----
b1_power$model <- factor(b1_power$model, levels = c("gecm", "adl", "adl22", "adl33"), labels =  c("gecm", "adl", "adl22", "adl33"))
ggplot(b1_power, aes(shape = model)) +
  geom_hline(yintercept = 1, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = b1_power), position = position_dodge(width = 0.75)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Power", x = "T", shape = "", title = expression("Power in the coefficient of x"[t])) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) +  
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("gecm", "adl", "adl22", "adl33"),
                     values = c(3, 4, 5, 6),
                     labels = c("GECM", "ADL(1,1)", "ADL(2,2)", "ADL(3,3)")) +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1) 
ggsave("GECM_b1_power.pdf")

#----Lx Bias----
b2_bias$model <- factor(b2_bias$model, levels = c("gecm", "adl", "adl22", "adl33"), labels =  c("gecm", "adl", "adl22", "adl33"))
ggplot(b2_bias, aes(shape = model)) +
  geom_hline(yintercept = 0, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = b2_bias), position = position_dodge(width = 0.75)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Bias", x = "T", shape = "", title = expression("Bias in the coefficient of x"[t-1])) + ylim(-0.5, 0.5) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("gecm", "adl", "adl22", "adl33"),
                     values = c(3, 4, 5, 6),
                     labels = c("GECM", "ADL(1,1)", "ADL(2,2)", "ADL(3,3)")) +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1) 
ggsave("GECM_b2_bias.pdf")

#----Lx Type 1----
b2_type1$model <- factor(b2_type1$model, levels = c("gecm", "adl", "adl22", "adl33"), labels =  c("gecm", "adl", "adl22", "adl33"))
ggplot(b2_type1, aes(shape = model)) +
  geom_hline(yintercept = 0.05, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = b2_type1), position = position_dodge(width = 0.75)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Type 1 Error Rate", x = "T", shape = "", title = expression("Type 1 error rate in the coefficient of x"[t-1])) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("gecm", "adl", "adl22", "adl33"),
                     values = c(3, 4, 5, 6),
                     labels = c("GECM", "ADL(1,1)", "ADL(2,2)", "ADL(3,3)")) +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1) 
ggsave("GECM_b2_type1.pdf")

#----Lx Coverage----
b2_type1 %>%
  mutate(b2_coverage = 1 - b2_type1) %>%
  ggplot(aes(shape = model)) +
  geom_hline(yintercept = 0.95, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = b2_coverage), position = position_dodge(width = 0.75)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Coverage Probability", x = "T", shape = "", title = expression("Coverage probability of the coefficient of x"[t-1])) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("gecm", "adl", "adl22", "adl33"),
                     values = c(3, 4, 5, 6),
                     labels = c("GECM", "ADL(1,1)", "ADL(2,2)", "ADL(3,3)")) +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1) 
ggsave("GECM_b2_coverage.pdf")

#----Lx Power----
b2_power$model <- factor(b2_power$model, levels = c("gecm", "adl", "adl22", "adl33"), labels =  c("gecm", "adl", "adl22", "adl33"))
ggplot(b2_power, aes(shape = model)) +
  geom_hline(yintercept = 1, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = b2_power), position = position_dodge(width = 0.75)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Power", x = "T", shape = "", title = expression("Power in the coefficient of x"[t-1])) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("gecm", "adl", "adl22", "adl33"),
                     values = c(3, 4, 5, 6),
                     labels = c("GECM", "ADL(1,1)", "ADL(2,2)", "ADL(3,3)")) +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1) 
ggsave("GECM_b2_power.pdf")

#----L2x Bias----
ggplot(b3_bias, aes(shape = model)) +
  geom_hline(yintercept = 0, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = b3_bias), position = position_dodge(width = 0.35)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Bias", x = "T", shape = "", title = expression("Bias in the coefficient of x"[t-2])) + ylim(-0.5, 0.5) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) +  
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("adl22", "adl33"),
                     values = c(5, 6),
                     labels = c("ADL(2,2)", "ADL(3,3)")) +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1) 
ggsave("GECM_b3_bias.pdf")

#----L2x Type 1----
ggplot(b3_type1, aes(shape = model)) +
  geom_hline(yintercept = 0.05, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = b3_type1), position = position_dodge(width = 0.35)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Type 1 Error Rate", x = "T", shape = "", title = expression("Type 1 error rate in the coefficient of x"[t-2])) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("adl22", "adl33"),
                     values = c(5, 6),
                     labels = c("ADL(2,2)", "ADL(3,3)")) +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1) 
ggsave("GECM_b3_type1.pdf")

#----L2x Coverage----
b3_type1 %>%
  mutate(b3_coverage = 1 - b3_type1) %>%
ggplot(aes(shape = model)) +
  geom_hline(yintercept = 0.95, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = b3_coverage), position = position_dodge(width = 0.35)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Coverage Probability", x = "T", shape = "", title = expression("Coverage probability of the coefficient of x"[t-2])) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("adl22", "adl33"),
                     values = c(5, 6),
                     labels = c("ADL(2,2)", "ADL(3,3)")) +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1) 
ggsave("GECM_b3_coverage.pdf")

#----L3x Bias----
ggplot(b4_bias, aes(shape = model)) +
  geom_hline(yintercept = 0, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = b4_bias), position = position_dodge(width = 0.35)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Bias", x = "T", shape = "", title = expression("Bias in the coefficient of x"[t-3])) + ylim(-0.5, 0.5) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10)) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("adl33"),
                     values = c(6),
                     labels = c("ADL(3,3)"))
ggsave("GECM_b4_bias.pdf")

#----L3x Type 1----
ggplot(b4_type1, aes(shape = model)) +
  geom_hline(yintercept = 0.05, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = b4_type1), position = position_dodge(width = 0.35)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Type 1 Error Rate", x = "T", shape = "", title = expression("Type 1 error rate in the coefficient of x"[t-3])) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10)) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("adl33"),
                     values = c(6),
                     labels = c("ADL(3,3)"))
ggsave("GECM_b4_type1.pdf")

#----L3x Coverage----
b4_type1 %>%
  mutate(b4_coverage = 1 - b4_type1) %>%
ggplot(aes(shape = model)) +
  geom_hline(yintercept = 0.95, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = b4_coverage), position = position_dodge(width = 0.35)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Coverage Probability", x = "T", shape = "", title = expression("Coverage probability of the coefficient of x"[t-3])) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10)) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("adl33"),
                     values = c(6),
                     labels = c("ADL(3,3)"))
ggsave("GECM_b4_coverage.pdf")

#----LRM Bias----
lrm_bias$model <- factor(lrm_bias$model, levels = c("gecm", "adl", "adl22", "adl33"), labels = c("gecm", "adl", "adl22", "adl33"))
ggplot(lrm_bias, aes(shape = model)) +
  geom_hline(yintercept = 0, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = lrm_bias), position = position_dodge(width = 0.75)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Bias", x = "T", shape = "", title = expression("Bias in the LRM of x"[t])) + ylim(-0.5, 0.5) + 
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) +  
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("gecm", "adl", "adl22", "adl33"),
                     values = c(3, 4, 5, 6),
                     labels = c("GECM", "ADL(1,1)", "ADL(2,2)", "ADL(3,3)")) +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1) 
ggsave("GECM_lrm_bias.pdf")

#----LRM Type 1----
lrm_type1$model <- factor(lrm_type1$model, levels = c("gecm", "adl", "adl22", "adl33"), labels = c("gecm", "adl", "adl22", "adl33"))
ggplot(lrm_type1, aes(shape = model)) +
  geom_hline(yintercept = 0.05, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = lrm_type1), position = position_dodge(width = 0.75)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Type 1 Error Rate", x = "T", shape = "", title = expression("Type 1 error rate in the LRM of x"[t])) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) +  
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("gecm", "adl", "adl22", "adl33"),
                     values = c(3, 4, 5, 6),
                     labels = c("GECM", "ADL(1,1)", "ADL(2,2)", "ADL(3,3)")) +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1) 
ggsave("GECM_lrm_type1.pdf")

#----LRM Coverage----
lrm_type1 %>%
  mutate(lrm_coverage = 1 - lrm_type1) %>%
  ggplot(aes(shape = model)) +
  geom_hline(yintercept = 0.95, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = lrm_coverage), position = position_dodge(width = 0.75)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Coverage Probability", x = "T", shape = "", title = expression("Coverage probability of the LRM of x"[t])) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) +  
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("gecm", "adl", "adl22", "adl33"),
                     values = c(3, 4, 5, 6),
                     labels = c("GECM", "ADL(1,1)", "ADL(2,2)", "ADL(3,3)")) +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1) 
ggsave("GECM_lrm_coverage.pdf")

#----LRM Power----
lrm_power$model <- factor(lrm_power$model, levels = c("gecm", "adl", "adl22", "adl33"), labels = c("gecm", "adl", "adl22", "adl33"))
ggplot(lrm_power, aes(shape = model)) +
  geom_hline(yintercept = 1, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = lrm_power), position = position_dodge(width = 0.75)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Power", x = "T", shape = "", title = expression("Power in the LRM of x"[t])) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("gecm", "adl", "adl22", "adl33"),
                     values = c(3, 4, 5, 6),
                     labels = c("GECM", "ADL(1,1)", "ADL(2,2)", "ADL(3,3)")) +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1) 
ggsave("GECM_lrm_power.pdf")

#----F-test, Lags of y Type 1 (sum = true phi_1)----
ftest_y_t1 %>%
  ggplot(aes(shape = model)) +
  geom_hline(yintercept = 0.05, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = ftest_y_t1), position = position_dodge(width = 0.35)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Type 1 Error Rate", x= "T", shape = "", title = expression("Type 1 error rate for a F-test on the coefficients of the lags of y"[t])) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("adl22", "adl33"),
                     values = c(5, 6),
                     labels = c("ADL(2,2)", "ADL(3,3)")) +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1)
ggsave("GECM_ftest_y_t1.pdf")

#----F-test, Lags of y Coverage (sum = true phi_1)----
ftest_y_t1 %>%
  mutate(ftest_y_coverage = 1 - ftest_y_t1) %>%
  ggplot(aes(shape = model)) +
  geom_hline(yintercept = 0.95, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = ftest_y_coverage), position = position_dodge(width = 0.35)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Coverage Probability", x= "T", shape = "", title = expression("Coverage probability for a F-test on the coefficients of the lags of y"[t])) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("adl22", "adl33"),
                     values = c(5, 6),
                     labels = c("ADL(2,2)", "ADL(3,3)")) +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1)
ggsave("GECM_ftest_y_coverage.pdf")

#----F-test, Lags of y Power (sum = 0)----
ftest_y_power %>%
  ggplot(aes(shape = model)) +
  geom_hline(yintercept = 1, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = ftest_y_power), position = position_dodge(width = 0.35)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Power", x= "T", shape = "", title = expression("Power for a F-test on the coefficients of the lags of y"[t])) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("adl22", "adl33"),
                     values = c(5, 6),
                     labels = c("ADL(2,2)", "ADL(3,3)")) +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1)
ggsave("GECM_ftest_y_power.pdf")

#----F-test, Lags of x Type 1 (sum = gamma)----
ftest_x_t1 %>%
  ggplot(aes(shape = model)) +
  geom_hline(yintercept = 0.05, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = ftest_x_t1), position = position_dodge(width = 0.35)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Type 1 Error Rate", x= "T", shape = "", title = expression("Type 1 error rate for a F-test on the coefficients of x"[t]~"and its lags")) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("adl22", "adl33"),
                     values = c(5, 6),
                     labels = c("ADL(2,2)", "ADL(3,3)")) +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1)
ggsave("GECM_ftest_x_t1.pdf")

#----F-test, Lags of x Coverage (sum = gamma)----
ftest_x_t1 %>%
  mutate(ftest_x_coverage = 1 - ftest_x_t1) %>%
  ggplot(aes(shape = model)) +
  geom_hline(yintercept = 0.95, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = ftest_x_coverage), position = position_dodge(width = 0.35)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Coverage Probability", x= "T", shape = "", title = expression("Coverage probability for a F-test on the coefficients of x"[t]~"and its lags")) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("adl22", "adl33"),
                     values = c(5, 6),
                     labels = c("ADL(2,2)", "ADL(3,3)")) +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1)
ggsave("GECM_ftest_x_coverage.pdf")

#----F-test, Lags of x Power (sum = 0)----
ftest_x_power %>%
  ggplot(aes(shape = model)) +
  geom_hline(yintercept = 1, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = ftest_x_power), position = position_dodge(width = 0.35)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Power", x= "T", shape = "", title = expression("Power for a F-test on the coefficients of x"[t]~"and its lags")) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("adl22", "adl33"),
                     values = c(5, 6),
                     labels = c("ADL(2,2)", "ADL(3,3)")) +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1)
ggsave("GECM_ftest_x_power.pdf")

